dat$InAgric_70<-dat$LForstwirtschaft_insg_1970/dat$Erwerbstaetige_insg_1970
dat$InManuf70 <-dat$Produzierenes_Gewerbe_insg_1970/dat$Erwerbstaetige_insg_1970
dat$HandelVerkehr70 <-dat$Handel_Verkehr_insg_1970/dat$Erwerbstaetige_insg_1970
dat$pAuslaender1970<-dat$Auslaender_insg_1970/dat$Bevoelkerung_insges_1970
dat$SelbstStandige_70<-dat$Selbstaendige_insg_1970/dat$Erwerbstaetige_insg_1970
dat$attendants_rel1970s<-dat$attendants_rel1970/100 #express percentage as share
dat$attendants_rel1980s<-dat$attendants_rel1980/100
dat$attendants_rel1990s<-dat$attendants_rel1990/100
dat$preligious_50s<-dat$preligious_50/100 #express percentage as share
dat$preligious_70s<-dat$preligious_70/100
dat$preligious_87s<-dat$preligious_87/100
dat$pdivorce_tomar_70<-dat$Geschieden_insg_1970/dat$Verheiratet_insg_1970
dat$pdivorce_tomar_87<-dat$Geschieden_insg_1987/dat$Verheiratet_insg_1987
#calculate party vote shares to eligible voters for CSU:
dat$csu_voters1949<-dat$votes_csu_1949/dat$eligible_voters_1949
dat$csu_voters1953<-dat$votes_csu_1953/dat$eligible_voters_1953
dat$csu_voters1957<-dat$votes_csu_1957/dat$eligible_voters_1957
dat$csu_voters1961<-dat$votes_csu_1961/dat$eligible_voters_1961
dat$csu_voters1965<-dat$votes_csu_1965/dat$eligible_voters_1965
dat$csu_voters1969<-dat$votes_csu_1969/dat$eligible_voters_1969
dat$csu_voters1972<-dat$votes_csu_1972/dat$eligible_voters_1972
dat$csu_voters1976<-dat$votes_csu_1976/dat$eligible_voters_1976
dat$csu_voters1980<-dat$votes_csu_1980/dat$eligible_voters_1980
dat$csu_voters1983<-dat$votes_csu_1983/dat$eligible_voters_1983
dat$csu_voters1987<-dat$votes_csu_1987/dat$eligible_voters_1987
dat$csu_voters1990<-dat$votes_csu_1990/dat$eligible_voters_1990
#calculate party vote shares to eligible voters for the Green Party:
dat$greens_voters1980<-dat$votes_greens_1980/dat$eligible_voters_1980
dat$greens_voters1983<-dat$votes_greens_1983/dat$eligible_voters_1983
dat$greens_voters1987<-dat$votes_greens_1987/dat$eligible_voters_1987
dat$greens_voters1990<-dat$votes_greens_1990/dat$eligible_voters_1990
#Create factor for fixed effects used in Conley error specifications
dat$county_id_f<-as.factor(dat$kreis_kenn1950_alt) #prewar kreise had same borders as 1950 kreise
dat$diocese_id_f<-as.factor(dat$diocese) #diocese fixed effect
###################################
######## FIGURE 2 in the paper ####
###################################
###########################################################################################################
# Replication for county-level analysis for Bavaria presented in Table A23 (sorting after resettlement)
# Charnysh and Schaub, Migration and Social Change: Evidence from post-WWII Displacement in Germany, JOP
# R version 4.3.1 (2023-06-16)
###########################################################################################################
rm(list = ls())
library(openxlsx)
library(stargazer)
library(sandwich)
library(lmtest)
dat3950<-read.xlsx("county_data/dataset19391950.xlsx") #file that includes distances and gender
dat3950$Evang1950[which(dat3950$KREIS50=="5536000")]<-13184  ##Correct Evang1950 for Landkreis Münster based on Gemeindestatistik Nordheim Westfalen, pages 1-23
mig6061<-read.xlsx("county_data/migration196061merged.xlsx")
dat<-merge(dat3950, mig6061, by="Kennziffer")
dat$pfemale1939<-dat$Women1939/dat$Tot1939
dat$ShareExp1950<-dat$Heimatvertriebene1950/dat$Pop1950
dat$ShareCath1939<-dat$Katholische1939/dat$Population1939
dat$ShareEvang1939<-dat$Evangelische1939/dat$Population1939
dat$ShareOther1939<-(dat$Population1939-dat$Evangelische1939-dat$Katholische1939)/dat$Population1939
dat$Frac1939<-1-(dat$ShareCath1939^2+dat$ShareEvang1939^2+dat$ShareOther1939^2)
summary(dat$Frac1939)
dat$ShareAgr1939<-dat$AgricultureInsgesamt1939/dat$Erwerbspersonen_Total1939
dat$ShareDestroyed<-dat$KriegsschaedenNWohnungen/dat$NormalWohnungen
dat$lnDistEastBorder<-log(dat$distance_to_borderkm)
dat$ShareCath1950<-dat$Kath1950/dat$Pop1950
dat[which(is.na(dat$KriegsschaedenNWohnungen)),] #Two places that are NAs in Braun's dataset.
dat$OtherReligion1950<-dat$Pop1950-dat$Kath1950-dat$Evang1950
summary(dat$OtherReligion1950)
dat$Frac1950<-1-((dat$Kath1950/dat$Pop1950)^2+(dat$Evang1950/dat$Pop1950)^2+(dat$OtherReligion1950/dat$Pop1950)^2)
summary(dat$Frac1950)
dat$lnPop1939<-log(dat$Population1939)
dat$PopDens1939<-dat$Population1939/dat$AREA50
dat$DeltaFrac<-dat$Frac1950-dat$Frac1939
#Calculate migration variables
dat$ptarrived1960<-dat$Zugezogene60/dat$Population1960
dat$ptleft1960<-dat$Fortgezogene60/dat$Population1960
dat$PtSurplusMig195061<-dat$SurplusInOutM19501961/dat$Pop1950 #Proportion of the original population that migrated in the next 10 years.
###########################################################################################################
# Replication for county-level analysis for Bavaria presented in Table A3 (sorting after resettlement)
# Charnysh and Schaub, Migration and Social Change: Evidence from post-WWII Displacement in Germany, JOP
# R version 4.3.1 (2023-06-16)
###########################################################################################################
rm(list = ls())
library(foreign)
library(openxlsx)
library(dplyr)
library(sandwich)
library(lmtest)
library(stargazer)
dat<-read.xlsx("county_data/dataset19391950.xlsx")
dat[which(dat$KREIS50=="5536000"), c("Pop1950", "Evang1950", "Kath1950", "Judische1950", "Heimatvertriebene1950", "VertrEvang1950", "VertrKath1950")]
#Correct Evang1950 for Landkreis Münster based on Gemeindestatistik Nordheim Westfalen, pages 1-23
dat$Evang1950[which(dat$KREIS50=="5536000")]<-13184
##Variables for 1950:
dat$srefugees<-dat$Heimatvertriebene1950/dat$Pop1950
dat$sprotexp<-dat$VertrEvang1950/dat$Heimatvertriebene1950 #Share Protestant expellees
dat$scathexp<-dat$VertrKath1950/dat$Heimatvertriebene1950 #Share Catholic expellees
dat$sotherexp<-(dat$Heimatvertriebene1950-dat$VertrEvang1950-dat$VertrKath1950)/dat$Heimatvertriebene1950
dat[which(dat$sotherexp<0), c("KREIS50", "Kreis","Heimatvertriebene1950", "VertrEvang1950", "VertrKath1950")] #some below zero, fix
dat$sotherexp<-ifelse(dat$sotherexp<0, 0, dat$sotherexp)  #Assumption that these are zeros.
dat$natives<-dat$Pop1950-dat$Heimatvertriebene1950 #
dat$OtherReligion1950<-dat$Pop1950-dat$Kath1950-dat$Evang1950
dat$scathnative<-(dat$Kath1950-dat$VertrKath1950)/dat$natives #Share Catholic natives
dat$sprotnative<-(dat$Evang1950-dat$VertrEvang1950)/dat$natives #Share Protestant natives
dat$sothernative<-(dat$OtherReligion1950-(dat$Heimatvertriebene1950-dat$VertrEvang1950-dat$VertrKath1950))/dat$natives #Share Other natives
dat[which(dat$sothernative<0),]  #all fine
dat$nativefrac1950<-1-(dat$scathnative^2+dat$sprotnative^2+dat$sothernative^2)
summary(dat$nativefrac1950)
dat$reffrac1950<-1-(dat$scathexp^2+dat$sprotexp^2+dat$sotherexp^2)
summary(dat$reffrac1950)
dat$sprotestants1939<-dat$Evangelische1939/dat$Population1939
dat$scatholics1939<-dat$Katholische1939/dat$Population1939
dat$OtherReligion1939<-dat$Population1939-dat$Katholische1939-dat$Evangelische1939
dat$sother1939<-dat$OtherReligion1939/dat$Population1939
dat$sprotestants1950<-dat$Evang1950/dat$Pop1950
dat$scatholics1950<-dat$Kath1950/dat$Pop1950
dat$sother1950<-(dat$Pop1950-dat$Kath1950-dat$Evang1950)/dat$Pop1950
dat$relfrac1939<-1-(dat$sprotestants1939^2+dat$scatholics1939^2+dat$sother1939^2)
summary(dat$relfrac1939)
dat$relfrac1950<-1-(dat$sprotestants1950^2+dat$scatholics1950^2+dat$sother1950^2)
dat$DeltaProt<-dat$sprotestants1950-dat$sprotestants1939
dat$DeltaCath<-dat$scatholics1950-dat$scatholics1939
dat$DeltaDiv<-dat$relfrac1950-dat$relfrac1939
##At the county level: did Catholic and Protestant shares of natives stay the same between 1939 and 1950? (reported on p. 12 in the paper)
cor.test(dat$relfrac1939, dat$nativefrac1950) #correlation of 0.98 (p<0.0001).
cor.test(dat$relfrac1939[dat$Land=="Bavaria"], dat$nativefrac1950[dat$Land=="Bavaria"]) #correlation of 0.98 (p<0.0001)
#calculate religious distance following Braun and Dwenger (2020):
dat$ReligDistanceV1<-sqrt((dat$scathnative-dat$scathexp)^2+(dat$sprotnative-dat$sprotexp)^2+(dat$sothernative-dat$sotherexp)^2) #+(data_1987m2$ShareOtherNatives-data_1987m2$ShareOtherExpellees)^2
summary(dat$ReligDistanceV1)
dat$ReligDistanceV1[which(dat$ReligDistanceV1>1)]<-1 #three are above one, cap them at one.
summary(dat$ReligDistanceV1)
##############################################################
######## Table A4: Undestanding denominational change  #######
##############################################################
m1<-lm(nativefrac1950~relfrac1939, data=dat)
m1.c <- coeftest(m1, vcov=vcovHC(m1, type="HC1"))
m2<-lm(nativefrac1950~relfrac1939+reffrac1950, data=dat)
m2.c <- coeftest(m2, vcov=vcovHC(m2, type="HC1"))
m3<-lm(nativefrac1950~relfrac1939+srefugees, data=dat)
m3.c <- coeftest(m3, vcov=vcovHC(m3, type="HC1"))
m4<-lm(relfrac1950~srefugees+nativefrac1950, data=dat)
m4.c <- coeftest(m4, vcov=vcovHC(m4, type="HC1"))
m5<-lm(relfrac1950~srefugees+ReligDistanceV1+nativefrac1950, data=dat)
m5.c <- coeftest(m5, vcov=vcovHC(m5, type="HC1"))
m6<-lm(relfrac1950~srefugees+ReligDistanceV1+I(ReligDistanceV1*srefugees)+nativefrac1950, data=dat)
m6.c <- coeftest(m6, vcov=vcovHC(m6, type="HC1"))
stargazer(m1, m2, m3, m4, m5, m6,
se = list(m1.c[, 2], m2.c[, 2], m3.c[, 2], m4.c[, 2], m5.c[, 2], m6.c[, 2]),
p = list(m1.c[, 4], m2.c[, 4], m3.c[, 4], m4.c[, 4], m5.c[, 4], m6.c[, 4]),
covariate.labels = c("Religious fractionalization (1939)",
"Expellee fractionalization (1950)",
"Share expellees (1950)",
"Religious distance (1950)",
"I(Religious distance* Share expellees)",
"Natives fractionalization (1950)"),
omit.stat = c("LL", "ser", "f"),
no.space=TRUE
)
###########################################################################################################
# Replication for correlations presented in Table A5 (Weimar electoral outcomes)
# Charnysh and Schaub, Migration and Social Change: Evidence from post-WWII Displacement in Germany, JOP
# R version 4.3.1 (2023-06-16)
###########################################################################################################
rm(list = ls())
library(openxlsx)
library(dplyr)
library(xtable)
dat_el<-read.xlsx("county_data/Weimar_elections.xlsx")
##Manually correct population for Landkreis Münster or use NA (typo in the book)
##Religious composition in 1939
dat_el$scatholics1939<-dat_el$Katholische1939/dat_el$Population1939
dat_el$sprotestants1939<-dat_el$Evangelische1939/dat_el$Population1939
dat_el$sother1939<-(dat_el$Population1939-dat_el$Evangelische1939-dat_el$Katholische1939)/dat_el$Population1939
dat_el$relfrac1939<-1-(dat_el$scatholics1939^2+dat_el$sprotestants1939^2+dat_el$sother1939^2)
##Religious composition in 1950:
dat_el$scatholics1950<-dat_el$Kath1950/dat_el$Pop1950
dat_el$sprotestants1950<-dat_el$Evang1950/dat_el$Pop1950
dat_el$sother1950<-(dat_el$Pop1950-dat_el$Evang1950-dat_el$Kath1950)/dat_el$Pop1950
dat_el$sother1950[which(dat_el$sother1950<0)]<-NA #Landkreis Münster. Total population is too small; there is a discrepancy in various statistical publications, replaced w/ NA
dat_el$sother1950[which(dat_el$sother1950<0)]<-NA #Landkreis Münster. Total population is too small; there is a discrepancy in various statistical publications, replaced w/ NA
dat_el$sother1950[which(dat_el$sother1950<0)]<-NA #Landkreis Münster. Total population is too small; there is a discrepancy in various statistical publications, replaced w/ NA
dat_el$relfrac1950<-1-(dat_el$scatholics1950^2+dat_el$sprotestants1950^2+dat_el$sother1950^2)
#Main explanatory variable:
dat_el$DeltaFrac<-dat_el$relfrac1950-dat_el$relfrac1939
## Electoral outcomes from the weimar period:
dat_el$SPD1928<-dat_el$n285spd/dat_el$n285wb
dat_el$DNVP1928<-dat_el$n285dnvp/dat_el$n285wb
dat_el$ZENTRUM1928<-dat_el$n285zx/dat_el$n285wb
dat_el$DDP1928<-dat_el$n285ddpx/dat_el$n285wb
dat_el$SPD1930<-dat_el$n309spd/dat_el$n309wb
dat_el$DNVP1930<-dat_el$n309dnvp/dat_el$n309wb
dat_el$ZENTRUM1930<-dat_el$n309zx/dat_el$n309wb
dat_el$DSTP1930<-dat_el$n309dstp/dat_el$n309wb
#July 1932 election:
dat_el$SPD1932j<-dat_el$n327spd/dat_el$n327wb
dat_el$DNVP1932j<-dat_el$n327dnvp/dat_el$n327wb
dat_el$ZENTRUM1932j<-dat_el$n327zx/dat_el$n327wb
dat_el$DSTP1932j<-dat_el$n327dstp/dat_el$n327wb
#November 1932 election:
dat_el$SPD1932n<-dat_el$n32nspd/dat_el$n32nwb
dat_el$DNVP1932n<-dat_el$n32ndnvp/dat_el$n32nwb
dat_el$ZENTRUM1932n<-dat_el$n32nzx/dat_el$n32nwb
dat_el$DSTP1932n<-dat_el$n32ndstp/dat_el$n327wb
# 1933 election:
dat_el$SPD1933<-dat_el$n333spd/dat_el$n333wb
dat_el$ZENTRUM1933 <-dat_el$n333zx/dat_el$n333wb
#Obtaining correlation coefficient and R2 between Delta diversity and electoral outcomes:
corrs<-rbind(cbind(cor.test(dat_el$DeltaFrac, dat_el$SPD1928)$estimate, summary(lm(DeltaFrac~ SPD1928, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$ZENTRUM1928)$estimate,summary(lm(DeltaFrac~ ZENTRUM1928, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$DNVP1928)$estimate,summary(lm(DeltaFrac~ DNVP1928, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$DDP1928)$estimate,summary(lm(DeltaFrac~ DDP1928, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$SPD1930)$estimate, summary(lm(DeltaFrac~ SPD1930, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$ZENTRUM1930)$estimate, summary(lm(DeltaFrac~ZENTRUM1930, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$DNVP1930)$estimate, summary(lm(DeltaFrac~DNVP1930, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$DSTP1930)$estimate, summary(lm(DeltaFrac~ DSTP1930, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$SPD1932j)$estimate, summary(lm(DeltaFrac~ SPD1932j, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$ZENTRUM1932j)$estimate, summary(lm(DeltaFrac~ZENTRUM1932j, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$DNVP1932j)$estimate, summary(lm(DeltaFrac~DNVP1932j, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$DSTP1932j)$estimate, summary(lm(DeltaFrac~ DSTP1932j, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$SPD1932n)$estimate, summary(lm(DeltaFrac~ SPD1932n, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$ZENTRUM1932n)$estimate, summary(lm(DeltaFrac~ZENTRUM1932n, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$DNVP1932n)$estimate, summary(lm(DeltaFrac~DNVP1932n, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$DSTP1932n)$estimate, summary(lm(DeltaFrac~ DSTP1932n, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$SPD1933)$estimate, summary(lm(DeltaFrac~ SPD1933, data=dat_el))$r.squared),
cbind(cor.test(dat_el$DeltaFrac, dat_el$ZENTRUM1933)$estimate, summary(lm(DeltaFrac~ZENTRUM1933, data=dat_el))$r.squared))
corrs<-as.data.frame(corrs, col.names=c("Variable","Pearson Correlation", "Proportion of explained variance (R^2)"),
row.names=c("SPD 1928",  "Zentrum 1928",  "DNVP 1928", "DDP 1928",
"SPD 1930","Zentrum 1930","DNVP 1930", "DSTP 1930",
"SPD 1932 (July)", "Zentrum 1932 (July)","DNVP 1932 (July)", "DSTP 1932 (July)",
"SPD 1932 (November)","Zentrum 1932 (November)",  "DNVP 1932 (November)", "DSTP 1932 (November)",
"SPD 1933", "Zentrum 1933"))
names(corrs)<-c("Pearson Correlation", "Propotion of explained variance (R^2)")
xtable(corrs, digits=3)
#####################################
# Replication for municipality-level analysis for
# Charnysh and Schaub, Migration and Social Change: Evidence from post-WWII Displacement in Germany, JOP
# R version 4.3.1 (2023-06-16)
#####################################
rm(list = ls())
# Load packages
load <- c("readstata13", "openxlsx", "sandwich", "lmtest", "stargazer", "dplyr", "conleyreg","sensemakr", "reshape2", "DirectEffects", "texreg","xtable", "interflex", "haven", "sf", "ggplot2", "geosphere", "reldist", "ape", "GWmodel")
lapply(load, library, character.only = TRUE)
####### CREATE FUNCTIONS ########
cse = function(reg) {
rob = sqrt(diag(vcovHC(reg, type = "HC1")))
return(rob)
}
#function to create formulas
formulas <- function(a, b) {
formulas <- list()
n <- 0
for (i in a) {
for (j in b) {
n <- n + 1
formulas[n] <- paste(i, j)
}
}
return(formulas)
}
#function to run models
run_models <- function(formlist) {
models <- list()
n <- 0
for (i in 1:length(formlist)){
n <- n + 1
models[[n]] <- eval(parse(text=formlist[[n]]))
}
return(models)
}
#vcov function
hc1<-function(x) vcovHC(x, type="HC1")
#function to extract coefficients + calculate errors
run_coef <- function(output) {
stdev <- list()
n <- 0
for (i in 1:length(output)){
n <- n + 1
stdev[[n]] <- coeftest(output[[n]], vcov. = hc1)
}
return(stdev)
}
#Function to format pvalues as stars for output from sequential g-estimation
pval_stars <- function(p) {
if (p < 0.01) {
return("***")
} else if (p < 0.05) {
return("**")
} else if (p < 0.1) {
return("*")
} else {
return("")
}
}
########################################
######## Compute main variables: #######
########################################
dat<-read.xlsx("municipal_data/municipality-dataset-short.xlsx")
###########################################################################################################
# Replication for county-level analysis presented in Table A13 (interdenominational marriages)
# Charnysh and Schaub, Migration and Social Change: Evidence from post-WWII Displacement in Germany, JOP
# R version 4.3.1 (2023-06-16)
###########################################################################################################
rm(list = ls())
# Load packages
load <- c("openxlsx", "sandwich", "lmtest", "stargazer")
lapply(load, library, character.only = TRUE)
#######################################################################################
## Table A.14. Religious diversity and the prevalence of Catholic-Protestant marriages
#######################################################################################
#Read in new data at stadtkreis level
dat<-read.xlsx("county_data/refugees_intermarriage_kreise.xlsx")
###########################################################################################################
# Replication for county-level analysis for Bavaria presented in Table A23 (sorting after resettlement)
# Charnysh and Schaub, Migration and Social Change: Evidence from post-WWII Displacement in Germany, JOP
# R version 4.3.1 (2023-06-16)
###########################################################################################################
rm(list = ls())
library(openxlsx)
library(stargazer)
library(sandwich)
library(lmtest)
dat3950<-read.xlsx("county_data/dataset19391950.xlsx") #file that includes distances and gender
dat3950$Evang1950[which(dat3950$KREIS50=="5536000")]<-13184  ##Correct Evang1950 for Landkreis Münster based on Gemeindestatistik Nordheim Westfalen, pages 1-23
mig6061<-read.xlsx("county_data/migration196061merged.xlsx")
dat<-merge(dat3950, mig6061, by="Kennziffer")
dat$pfemale1939<-dat$Women1939/dat$Tot1939
dat$ShareExp1950<-dat$Heimatvertriebene1950/dat$Pop1950
dat$ShareCath1939<-dat$Katholische1939/dat$Population1939
dat$ShareEvang1939<-dat$Evangelische1939/dat$Population1939
dat$ShareOther1939<-(dat$Population1939-dat$Evangelische1939-dat$Katholische1939)/dat$Population1939
dat$Frac1939<-1-(dat$ShareCath1939^2+dat$ShareEvang1939^2+dat$ShareOther1939^2)
summary(dat$Frac1939)
dat$ShareAgr1939<-dat$AgricultureInsgesamt1939/dat$Erwerbspersonen_Total1939
dat$ShareDestroyed<-dat$KriegsschaedenNWohnungen/dat$NormalWohnungen
dat$lnDistEastBorder<-log(dat$distance_to_borderkm)
dat$ShareCath1950<-dat$Kath1950/dat$Pop1950
dat[which(is.na(dat$KriegsschaedenNWohnungen)),] #Two places that are NAs in Braun's dataset.
dat$OtherReligion1950<-dat$Pop1950-dat$Kath1950-dat$Evang1950
summary(dat$OtherReligion1950)
dat$Frac1950<-1-((dat$Kath1950/dat$Pop1950)^2+(dat$Evang1950/dat$Pop1950)^2+(dat$OtherReligion1950/dat$Pop1950)^2)
summary(dat$Frac1950)
dat$lnPop1939<-log(dat$Population1939)
dat$PopDens1939<-dat$Population1939/dat$AREA50
dat$DeltaFrac<-dat$Frac1950-dat$Frac1939
#Calculate migration variables
dat$ptarrived1960<-dat$Zugezogene60/dat$Population1960
dat$ptleft1960<-dat$Fortgezogene60/dat$Population1960
dat$PtSurplusMig195061<-dat$SurplusInOutM19501961/dat$Pop1950 #Proportion of the original population that migrated in the next 10 years.
################################################################################
############# Table A24: Religious diversity and in/out migration  #############
################################################################################
reg1<-lm(PtSurplusMig195061~DeltaFrac+Frac1939, data=dat)
reg2<-lm(PtSurplusMig195061~DeltaFrac+Frac1939+ShareExp1950+pfemale1939+ShareAgr1939+ShareDestroyed+lnDistEastBorder+PopDens1939+lnPop1939, data=dat)
reg3<-lm(ptarrived1960~DeltaFrac+Frac1939, data=dat)
reg4<-lm(ptarrived1960~DeltaFrac +Frac1939+ShareExp1950+pfemale1939+ShareAgr1939+ShareDestroyed+lnDistEastBorder+PopDens1939+lnPop1939, data=dat)
reg5<-lm(ptleft1960~DeltaFrac+Frac1939, data=dat)
reg6<-lm(ptleft1960~DeltaFrac +Frac1939+ShareExp1950+pfemale1939+ShareAgr1939+ShareDestroyed+lnDistEastBorder+PopDens1939+lnPop1939, data=dat)
ses1 <- list(coeftest(reg1, vcov = vcovHC(reg1, type="HC1"))[,2], coeftest(reg2, vcov = vcovHC(reg2, type="HC1"))[,2],
coeftest(reg3, vcov = vcovHC(reg3, type="HC1"))[,2], coeftest(reg4, vcov = vcovHC(reg4, type="HC1"))[,2],
coeftest(reg5, vcov = vcovHC(reg5, type="HC1"))[,2], coeftest(reg6, vcov = vcovHC(reg6, type="HC1"))[,2])
pvals1 <- list(coeftest(reg1, vcov = vcovHC(reg1, type="HC1"))[,4], coeftest(reg2, vcov = vcovHC(reg2, type="HC1"))[,4],
coeftest(reg3, vcov = vcovHC(reg3, type="HC1"))[,4], coeftest(reg4, vcov = vcovHC(reg4, type="HC1"))[,4],
coeftest(reg5, vcov = vcovHC(reg5, type="HC1"))[,4], coeftest(reg6, vcov = vcovHC(reg6, type="HC1"))[,4])
means_dv <- c(rep(round(mean(dat$PtSurplusMig195061,na.rm=T),2), 2),rep(round(mean(dat$ptarrived1960,na.rm=T),2), 2), rep(round(mean(dat$ptleft1960,na.rm=T),2), 2))
sds_dv <- c(rep(round(sd(dat$PtSurplusMig195061,na.rm=T),2), 2), rep(round(sd(dat$ptarrived1960,na.rm=T),2), 2),rep(round(sd(dat$ptleft1960,na.rm=T),2), 2))
mean_row <- c("DV mean", means_dv)
sd_row <- c("DV sd", sds_dv)
stargazer(
reg1, reg2, reg3, reg4, reg5, reg6,
se = ses1, p = pvals1, digits = 2,
covariate.labels = c(
"Delta diversity", "Fractionalization (1939)", "Share expellees",
"Share female (1939)", "Share in agriculture (1939)", "Share destroyed",
"Ln dist to the Eastern border", "Pop density (1939)", "Ln population (1939)"
),
omit = c("Constant"), omit.stat = c("ser", "f", "rsq"), no.space = TRUE,
column.labels = c("Net in/out migration 1950-1961", "In-migration 1960", "Out-migration 1960"),
title = "Net in/out migration. Robust SEs in parentheses.",
add.lines = list(mean_row, sd_row)  # Add the custom rows
)
###########################################################################################################
# Replication for county-level analysis for Bavaria presented in Table A24 (sorting after resettlement)
# Charnysh and Schaub, Migration and Social Change: Evidence from post-WWII Displacement in Germany, JOP
# R version 4.3.1 (2023-06-16)
###########################################################################################################
rm(list = ls())
library(foreign)
library(openxlsx)
library(dplyr)
library(sandwich)
library(lmtest)
library(stargazer)
dat<-read.xlsx("county_data/dataset19391950.xlsx")
dat[which(dat$KREIS50=="5536000"), c("Pop1950", "Evang1950", "Kath1950", "Judische1950", "Heimatvertriebene1950", "VertrEvang1950", "VertrKath1950")]
#Correct Evang1950 for Landkreis Münster based on Gemeindestatistik Nordheim Westfalen, pages 1-23
dat$Evang1950[which(dat$KREIS50=="5536000")]<-13184
##Variables for 1950:
dat$srefugees<-dat$Heimatvertriebene1950/dat$Pop1950
dat$sprotexp<-dat$VertrEvang1950/dat$Heimatvertriebene1950 #Share Protestant expellees
dat$scathexp<-dat$VertrKath1950/dat$Heimatvertriebene1950 #Share Catholic expellees
dat$sotherexp<-(dat$Heimatvertriebene1950-dat$VertrEvang1950-dat$VertrKath1950)/dat$Heimatvertriebene1950
dat[which(dat$sotherexp<0), c("KREIS50", "Kreis","Heimatvertriebene1950", "VertrEvang1950", "VertrKath1950")] #some below zero, fix
dat$sotherexp<-ifelse(dat$sotherexp<0, 0, dat$sotherexp)  #Assumption that these are zeros.
dat$natives<-dat$Pop1950-dat$Heimatvertriebene1950 #
dat$OtherReligion1950<-dat$Pop1950-dat$Kath1950-dat$Evang1950
dat$scathnative<-(dat$Kath1950-dat$VertrKath1950)/dat$natives #Share Catholic natives
dat$sprotnative<-(dat$Evang1950-dat$VertrEvang1950)/dat$natives #Share Protestant natives
dat$sothernative<-(dat$OtherReligion1950-(dat$Heimatvertriebene1950-dat$VertrEvang1950-dat$VertrKath1950))/dat$natives #Share Other natives
dat[which(dat$sothernative<0),]  #all fine
dat$nativefrac1950<-1-(dat$scathnative^2+dat$sprotnative^2+dat$sothernative^2)
summary(dat$nativefrac1950)
dat$reffrac1950<-1-(dat$scathexp^2+dat$sprotexp^2+dat$sotherexp^2)
summary(dat$reffrac1950)
dat$sprotestants1939<-dat$Evangelische1939/dat$Population1939
dat$scatholics1939<-dat$Katholische1939/dat$Population1939
dat$OtherReligion1939<-dat$Population1939-dat$Katholische1939-dat$Evangelische1939
dat$sother1939<-dat$OtherReligion1939/dat$Population1939
dat$sprotestants1950<-dat$Evang1950/dat$Pop1950
dat$scatholics1950<-dat$Kath1950/dat$Pop1950
dat$sother1950<-(dat$Pop1950-dat$Kath1950-dat$Evang1950)/dat$Pop1950
dat$relfrac1939<-1-(dat$sprotestants1939^2+dat$scatholics1939^2+dat$sother1939^2)
summary(dat$relfrac1939)
dat$relfrac1950<-1-(dat$sprotestants1950^2+dat$scatholics1950^2+dat$sother1950^2)
dat$DeltaProt<-dat$sprotestants1950-dat$sprotestants1939
dat$DeltaCath<-dat$scatholics1950-dat$scatholics1939
dat$DeltaDiv<-dat$relfrac1950-dat$relfrac1939
##At the county level: did Catholic and Protestant shares of natives stay the same between 1939 and 1950? (reported on p. 12 in the paper)
cor.test(dat$relfrac1939, dat$nativefrac1950) #correlation of 0.98 (p<0.0001).
cor.test(dat$relfrac1939[dat$Land=="Bavaria"], dat$nativefrac1950[dat$Land=="Bavaria"]) #correlation of 0.98 (p<0.0001)
#calculate religious distance following Braun and Dwenger (2020):
dat$ReligDistanceV1<-sqrt((dat$scathnative-dat$scathexp)^2+(dat$sprotnative-dat$sprotexp)^2+(dat$sothernative-dat$sotherexp)^2) #+(data_1987m2$ShareOtherNatives-data_1987m2$ShareOtherExpellees)^2
summary(dat$ReligDistanceV1)
dat$ReligDistanceV1[which(dat$ReligDistanceV1>1)]<-1 #three are above one, cap them at one.
summary(dat$ReligDistanceV1)
###########################################################################################################
# Replication for county-level analysis for Bavaria presented in Table A4 (Religious diversity in 1939 and 1950)
# Charnysh and Schaub, Migration and Social Change: Evidence from post-WWII Displacement in Germany, JOP
# R version 4.3.1 (2023-06-16)
###########################################################################################################
rm(list = ls())
library(foreign)
library(openxlsx)
library(dplyr)
library(sandwich)
library(lmtest)
library(stargazer)
dat<-read.xlsx("county_data/dataset19391950.xlsx")
dat[which(dat$KREIS50=="5536000"), c("Pop1950", "Evang1950", "Kath1950", "Judische1950", "Heimatvertriebene1950", "VertrEvang1950", "VertrKath1950")]
#Correct Evang1950 for Landkreis Münster based on Gemeindestatistik Nordheim Westfalen, pages 1-23
dat$Evang1950[which(dat$KREIS50=="5536000")]<-13184
##Variables for 1950:
dat$srefugees<-dat$Heimatvertriebene1950/dat$Pop1950
dat$sprotexp<-dat$VertrEvang1950/dat$Heimatvertriebene1950 #Share Protestant expellees
dat$scathexp<-dat$VertrKath1950/dat$Heimatvertriebene1950 #Share Catholic expellees
dat$sotherexp<-(dat$Heimatvertriebene1950-dat$VertrEvang1950-dat$VertrKath1950)/dat$Heimatvertriebene1950
dat[which(dat$sotherexp<0), c("KREIS50", "Kreis","Heimatvertriebene1950", "VertrEvang1950", "VertrKath1950")] #some below zero, fix
dat$sotherexp<-ifelse(dat$sotherexp<0, 0, dat$sotherexp)  #Assumption that these are zeros.
dat$natives<-dat$Pop1950-dat$Heimatvertriebene1950 #
dat$OtherReligion1950<-dat$Pop1950-dat$Kath1950-dat$Evang1950
dat$scathnative<-(dat$Kath1950-dat$VertrKath1950)/dat$natives #Share Catholic natives
dat$sprotnative<-(dat$Evang1950-dat$VertrEvang1950)/dat$natives #Share Protestant natives
dat$sothernative<-(dat$OtherReligion1950-(dat$Heimatvertriebene1950-dat$VertrEvang1950-dat$VertrKath1950))/dat$natives #Share Other natives
dat[which(dat$sothernative<0),]  #all fine
dat$nativefrac1950<-1-(dat$scathnative^2+dat$sprotnative^2+dat$sothernative^2)
summary(dat$nativefrac1950)
dat$reffrac1950<-1-(dat$scathexp^2+dat$sprotexp^2+dat$sotherexp^2)
summary(dat$reffrac1950)
dat$sprotestants1939<-dat$Evangelische1939/dat$Population1939
dat$scatholics1939<-dat$Katholische1939/dat$Population1939
dat$OtherReligion1939<-dat$Population1939-dat$Katholische1939-dat$Evangelische1939
dat$sother1939<-dat$OtherReligion1939/dat$Population1939
dat$sprotestants1950<-dat$Evang1950/dat$Pop1950
dat$scatholics1950<-dat$Kath1950/dat$Pop1950
dat$sother1950<-(dat$Pop1950-dat$Kath1950-dat$Evang1950)/dat$Pop1950
dat$relfrac1939<-1-(dat$sprotestants1939^2+dat$scatholics1939^2+dat$sother1939^2)
summary(dat$relfrac1939)
dat$relfrac1950<-1-(dat$sprotestants1950^2+dat$scatholics1950^2+dat$sother1950^2)
dat$DeltaProt<-dat$sprotestants1950-dat$sprotestants1939
dat$DeltaCath<-dat$scatholics1950-dat$scatholics1939
dat$DeltaDiv<-dat$relfrac1950-dat$relfrac1939
##At the county level: did Catholic and Protestant shares of natives stay the same between 1939 and 1950? (reported on p. 12 in the paper)
cor.test(dat$relfrac1939, dat$nativefrac1950) #correlation of 0.98 (p<0.0001).
cor.test(dat$relfrac1939[dat$Land=="Bavaria"], dat$nativefrac1950[dat$Land=="Bavaria"]) #correlation of 0.98 (p<0.0001)
#calculate religious distance following Braun and Dwenger (2020):
dat$ReligDistanceV1<-sqrt((dat$scathnative-dat$scathexp)^2+(dat$sprotnative-dat$sprotexp)^2+(dat$sothernative-dat$sotherexp)^2) #+(data_1987m2$ShareOtherNatives-data_1987m2$ShareOtherExpellees)^2
summary(dat$ReligDistanceV1)
dat$ReligDistanceV1[which(dat$ReligDistanceV1>1)]<-1 #three are above one, cap them at one.
summary(dat$ReligDistanceV1)
##############################################################
######## Table A4: Undestanding denominational change  #######
##############################################################
m1<-lm(nativefrac1950~relfrac1939, data=dat)
m1.c <- coeftest(m1, vcov=vcovHC(m1, type="HC1"))
m2<-lm(nativefrac1950~relfrac1939+reffrac1950, data=dat)
m2.c <- coeftest(m2, vcov=vcovHC(m2, type="HC1"))
m3<-lm(nativefrac1950~relfrac1939+srefugees, data=dat)
m3.c <- coeftest(m3, vcov=vcovHC(m3, type="HC1"))
m4<-lm(relfrac1950~srefugees+nativefrac1950, data=dat)
m4.c <- coeftest(m4, vcov=vcovHC(m4, type="HC1"))
m5<-lm(relfrac1950~srefugees+ReligDistanceV1+nativefrac1950, data=dat)
m5.c <- coeftest(m5, vcov=vcovHC(m5, type="HC1"))
m6<-lm(relfrac1950~srefugees+ReligDistanceV1+I(ReligDistanceV1*srefugees)+nativefrac1950, data=dat)
m6.c <- coeftest(m6, vcov=vcovHC(m6, type="HC1"))
stargazer(m1, m2, m3, m4, m5, m6,
se = list(m1.c[, 2], m2.c[, 2], m3.c[, 2], m4.c[, 2], m5.c[, 2], m6.c[, 2]),
p = list(m1.c[, 4], m2.c[, 4], m3.c[, 4], m4.c[, 4], m5.c[, 4], m6.c[, 4]),
covariate.labels = c("Religious fractionalization (1939)",
"Expellee fractionalization (1950)",
"Share expellees (1950)",
"Religious distance (1950)",
"I(Religious distance* Share expellees)",
"Natives fractionalization (1950)"),
omit.stat = c("LL", "ser", "f"),
no.space=TRUE
)
